/*********************************************************************************
 * Software License Agreement (BSD License)
 * Copyright (c) 2012, Carlos D. Correa (info@ngraph.org)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of netzen nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************************/
#include "timing.h"

#include "node.h"
#include "centrality.h"

#include <math.h>
#include <utility>                   // for std::pair
#include <algorithm>                 // for std::for_each
#include <limits.h>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/johnson_all_pairs_shortest.hpp>
//#include <boost/vector_property_map.hpp>

using namespace boost;
  
Closeness::Closeness(Controller *controller):AnalysisTool("Centrality","Closeness",controller) {
    setPropertyBool("Derivatives", false);
    setPropertyStringList("Edge Weight", "Default", "edge");
    setPropertyString("Result", this->getName());
}

void Closeness::compute() {
  assert(controller);
  // TO DO: set up these two via UI
  Graph *g = controller->getGraph("graph");
  bool derivs =  this->getPropertyBool("Derivatives");

  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";

  compute(g, derivs);
  controller->addDataProperty("node", resultProp, PROPERTY_FLOAT);
  g->computeNodeMinMaxFloat(resultProp);
  if(derivs) {
    controller->addDataProperty("edge",resultDerivs, PROPERTY_FLOAT);
    g->computeEdgeMinMaxFloat(resultDerivs);
  }
}

 
int Closeness::compute(Graph *graph, bool needDerivs) {
  timestamp t1 = now();
  int n = (int) graph->numNodes();
  
  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";
  string edgeProperty = getPropertyString("Edge Weight");
  if(listener) listener->setProgress(0);

  //  typedef adjacency_list<vecS, vecS, directedS> Graph;
  typedef adjacency_list<vecS, vecS, directedS, no_property,
      property< edge_weight_t, float, 
      property< edge_weight2_t, float > > > Graph;
   Graph g;

  int numEdges=0;
  vector<float> weights;
	ElementIterator *it = graph->createEdgeIterator();
	for(it->begin();it->hasMoreElements();) {
		Edge *edge = (Edge*) it->next();
		add_edge(edge->getNode0()->getIndex(), edge->getNode1()->getIndex(), g);
		float weight = graph->getAsFloat(edge, edgeProperty);
		weights.push_back(weight);
		numEdges++;
  }
	graph->destroyIterator(it);

  property_map < Graph, edge_weight_t >::type w = get(edge_weight, g);
  float *weights0 = new float[numEdges];
  float *wp = weights0;
  for(int i=0;i<numEdges;i++) wp[i] = weights[i];

  graph_traits < Graph >::edge_iterator e, e_end;
  for (boost::tie(e, e_end) = edges(g); e != e_end; ++e) {
    w[*e] = *wp++;
  }
  
  std::vector < float >d(n, (std::numeric_limits < float >::max)());
  float **D;
  D = new float*[n];
  for(int i=0;i<n;i++) D[i] = new float[n];
  johnson_all_pairs_shortest_paths(g, D, distance_map(&d[0]));


  vector<float> closeness;
  for(int i=0;i<n;i++) {
    float sum = 0;
    for(int j=0;j<n;j++) {
      sum+=D[i][j];
    }
    float close = (fabs(sum)>0)? 1.0/sum: 0.0;
    graph->getNode(i)->setPropertyFloat(resultProp, close);
    closeness.push_back(close);
  }
  

  timestamp t2 = now();
  if(needDerivs) {
    std::vector < float >d2(n, (std::numeric_limits < float >::max)());
    float **D2;
    D2 = new float*[n];
    for(int i=0;i<n;i++) {
      D2[i] = new float[n];
      for(int j=0;j<n;j++) {
				D2[i][j] = 0;
      }
    }
		double *sumWeights = new double[n];
		for(int gi=0;gi<n;gi++) {
			float sum = 0;
			ElementIterator *it = graph->getEdges(gi);
			for(it->begin();it->hasMoreElements();) {
				Edge *edge = (Edge*) it->next();
				float weight = graph->getAsFloat(edge, edgeProperty);
				sum += weight;
			}
			graph->destroyIterator(it);
			sumWeights[gi] = sum;
		}
    
    for(int i=0;i<n;i++) {
      //print("Computing Derivative w.r.t %d\n", i);
      if(listener) listener->setProgress((float) i/(float) n);
      int numEdges = 0;
			ElementIterator *it = graph->createEdgeIterator();
			for(it->begin();it->hasMoreElements();) {
				Edge *edge = (Edge*) it->next();
				float weight = graph->getAsFloat(edge,edgeProperty);
				int src = edge->getNode0()->getIndex();
				int dst = edge->getNode1()->getIndex();
				weights[numEdges] = ((i==src || i==dst) && sumWeights[src]!=0)? weight + weight/sumWeights[src]: weight;
				numEdges++;
      }
      graph_traits < Graph >::edge_iterator e, e_end;
      int k =0;
      for (boost::tie(e, e_end) = edges(g); e != e_end; ++e) {
				w[*e] = weights[k++];
      }
      
      johnson_all_pairs_shortest_paths(g, D2, distance_map(&d2[0]));
      for(int j=0;j<n;j++) {
				float sum = 0;
				for(int k=0;k<n;k++) {
					sum+=D2[j][k];
				}
				float close = (sum>0)? 1.0/sum: 0.0;
				bool transpose = true;
				if(transpose) {
					if(graph->hasEdge(i,j)) {
						graph->getEdge(i,j)->setPropertySigned(resultDerivs, close-closeness[j]);
					}
				} else {
					if(graph->hasEdge(j,i)) {
						graph->getEdge(j,i)->setPropertySigned(resultDerivs, close-closeness[j]);
					}
				}
      }
      //print("Computed derivative w.r.t %d\n", i);
    }
		delete sumWeights;
    delete[] D2;
  }

  timestamp t3 = now();
  if(listener) listener->done();
 
  delete[] D;
  printf("[CLOSENESS] Time  %f sec\n", t2 - t1);
  printf("[CLOSENESS DERIVS] Time  %f sec\n", t3 - t2);
  delete weights0;

  return 1;
}


